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Abstract 

We derive the exact equations of motion for the random neighbor version 
of the Olami-Feder-Christensen earthquake model in the infinite-size limit. We 
solve them numerically, and compare with simulations of the model for large 
numbers of sites. We find perfect agreement. But we do not find any scaling 
or phase transitions, except in the conservative limit. This is in contradiction 
to claims by Lise & Jensen (Phys. Rev. Lett. 76, 2326 (1996)) based on 
approximate solutions of the same model. It indicates again that scaling in the 
Olami-Feder-Christensen model is only due to partial synchronization driven by 
spatial inhomogeneities. Finally, we point out that our method can be used also 
for other SOC models, and treat in detail the random neighbor version of the 
Feder-Feder model. 
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1 Introduction 



During the last ten years more than 2000 pubhcations were concerned with the idea 
of self-organized criticality (SOC) proposed by Bak, Tang and Wiesenfeld(BTW) 
They introduced a non-equilibrium system, the so-called sandpile model, which is 
driven slowly by adding single sand grains at random positions. Without any control 
parameter to fine-tune, it evolves into a critical state. In this state the system reacts 
to the external drive with a series of relaxation events (avalanches) . It becomes critical 
in the sense that the spatial and temporal distributions of these avalanches obey power 
laws, indicating that any characteristic scales in space and time are lost. The attribute 
'self-organized' is to stress the absence of a fine-tuned control parameter. 

A crucial point in understanding the robust scaling of the BTW model is the exis- 
tence of a conservation law ^j: the total amount of sand in the system is conserved, 
if boundary effects and external perturbations are neglected. 

In the frame of this concept, Olami, Feder, and Christensen (OFC) introduced 
a nonconservative 'continuous cellular automaton' [0] as a specific realization of the 
two-dimensional Burridge-Knopoff earthquake model 0]. Details of this model will be 
described below. In contrast to the BTW model it is not conservative in general. It 
involves a parameter a, and a conservation law holds only for a specific value a = ac- 
It was found in [0, ^ and subsequent simulations that the system displays power 
law behavior in a wide range of the control parameter a (not only near ac), and the 
critical exponents depend on a. Thus the model seems to show SOC, and conservation 
seems not a necessary condition. 

But, on the other hand, it seems that spatial inhomogeneities are crucial for the 
observation of scaling in the OFC model [§], 0; 3- In the original paper by OFC the 
boundary conditions (be) were not periodic, which induced an inhomogeneity with 
a diverging length scale in the thermodynamic limit. This inhomogeneity of the be 
leads to partial synchronization in the bulk which is both driven and destroyed by the 
boundary [Q. Subsequent simulations with periodic be showed no scaling 0, 0, as did 



also simulations with frozen randomness but without diverging length scales |jl0|, |Tl 
The basic source of scaling in the OFC model is the slow build-up of large coherent 
domains in which the system itself is homogeneous, but which are driven by regions 
where the system is not homogeneous. 

Although the definitions of these SOC models are simple, and they are easily sim- 
ulated on a computer, only few exact results are known. Most of the difficulties in 
the analytical treatment arise from the spatial correlations due to the interactions of 
the particles. In a mean-field theory, which is the first step towards a detailed un- 
derstanding, these correlations are simply neglected. A more refined strategy to avoid 
spatial correlations is to replace the nearest neighbors interactions by interactions be- 
tween random sites. For the OFC model this was already attempted by |jl2|. But in 



that paper additional assumptions and approximations were made which are hard to 
justify. With these assumptions, a transition was found from non-SOC to SOC at 
a significantly less than ac- This is very surprising, as we argued above that spatial 
structures are crucial for the emergence of scaling, and any such structures are of course 
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eliminated in the random neighbor version. 

In the following we study the random neighbor model in detail without any further 
approximations. We will be led to a complete set of equations which allow us to 
calculate numerically all the relevant quantities. We will see that there is no SOC 
in the dissipative regime of the control-parameter. In the case of conservation the 
exact solution shows that the system becomes a critical branching process equivalent 
to critical percolation on a Bethe tree, and the critical exponents take their mean-field 
values. 

2 The model 

The model lives on a set of sites, each of them equipped with a continuous stress 
(or 'force') variable Zj. Each Zi can take any value > 0, but only values < 1 are stable. 
After having initialized each site with a randomly chosen value Zi G [0, 1[, the system 
evolves according to the following rules: 

(i) All Zi are simultaneously and continuously increased with the same speed v = 1. 

(ii) If any Zi exceeds the threshold value 1 the above driving stops, and the forces 
are redistributed in the following way: 

All unstable sites discharge simultaneously, 

^. _^ Vz, > 1 . (1) 

For each of these discharging sites n random "neighbors", are chosen and 

their stress variables are increased by a fixed fraction of Zi. 

Zj^ — >• Zji^ + azi k = 1, . . . ,n (2) 

The integer n is constant but otherwise arbitrary. If the application of eq. (0) creates 
new unstable sites, rule (ii) is again applied in the next time step, again simultaneously 
for all unstable sites. This procedure is repeated until all sites are stable. After that, 
the system is again driven according to rule (i), until at least one site with Zi = 1 
appears. A series of causally connected discharging events is called an earthquake 
or an avalanche. Its size s is measured by the total number of discharges. If a site 
discharges m times during an avalanche, it is counted m-fold in the calculation of s. 
The duration t of an earthquake is defined as the number of sweeps through the lattice 
necessary to get a stable configuration. Obviously s as well as t is always > 1. 

The parameter a which controls the dissipation can take any value between and 
l/n(Q!>l/nis unphysical since sooner or later an infinite and ever growing avalanche 
would occur). Only for a = 1/n the system is conservative. Note that the randomness 
of the neighbors appearing in eq. (0) is annealed: For each discharging event, the n ran- 
dom neighbors are chosen anew. Obviously this prevents that any spatial correlations 
in the values of z to build up. 

The numerical calculations as well as the simulations are restricted to the case n = 4. 
Obviously this is most appropriate for a mean-field theory of the two-dimensional OFC- 
model. But our analytic results are more general and hold for any n >2. 
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As well known, the original (nearest neighbor) version of the model is very sensitive 
to the choice of bc's 0. Any be other than periodic introduce inhomogeneities 
which are crucial in building up the spatial structures which manifest themselves in 
non-trivial avalanches |^ Q . For the random neighbor version, non-periodic be were 
used in [|^. This also introduces spatial inhomogeneity which is however completely 
irrelevant for the dynamics, 'space' being a dummy concept in a random neighbor 



model. In addition, the be used in [|12| lead to specific finite size corrections which 



might be not easy to disentangle from the true asymptotic behavior. In contrast, we 
treat all sites equally in the present paper, mimicking thereby periodic be. In addition, 
we shall only study the infinite size limit. More precisely, we shall formally work with 
a finite number N of sites, but will understand that we are only interested in the limit 
^ oo. For finite sizes there are correlations which make the study of the model 
rather awkward. 



3 Random Neighbor Theory 

In the OFC model, there is a finite chance that two sites become unstable simultane- 
ously during the continuous increase (i). It arises from the non-zero probability that 
two sites which had discharged in the same previous earth quake have not been hit by 
a discharging neighbor (or have been hit by the same neighbors) until they reach z = 1. 
In the lattice version this implies that the notion of an earth quake itself becomes a bit 
delicate: should we consider an event which was triggered simultaneously by two sites 
as one earth quake or two? In the off-lattice version we still have a non-zero chance 
for such events. But on an infinite lattice the sub-quakes following each unstable site 
will not overlap. Thus they will evolve completely independently. This means that the 
model becomes effectively abelian |]13| in the sense that we can change the order of 
updates in different sub-events. Also, we can associate earth quakes uniquely with the 
original unstable sites which triggered them. In the following, we will always define 
earth quakes in this way. An event which started with k sites becoming unstable is 
counted as k earth quakes, separated by infinitesimal time delays and taking place in 
arbitrary order. 

So we can assume without loss of generality that after the relaxation of an earth- 
quake there is exactly one site which has a stress value greater than all the others, 
and which will be the seed of the next avalanche. The value of this stress immediately 
after the earthquake has stopped will be called z^. Its mean value, averaged over all 
earthquakes, is denoted by z^. Since we consider the large system limit, z^ will not be 
correlated with the size of the previous avalanche. This is our crucial assumption, and 
it depends on the fact that we can neglect 'global' avalanches whose size is comparable 
to the total size of the system. In this limit the model thus becomes a branching pro- 
cess with time dependent branching rates. We shall later verify that this assumption 
is self consistent, and is true in simulations. 

The average increase of the force on each of the N sites between two earth quakes 
due to the external driving is then given by 1 — 'z^. On the other hand, each discharge 
dissipates an average value of {l—na)z^^rist, where z^^st is the mean force on the unstable 
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sites, averaged over all discharging events. In the stationary state, when Zi fluctuates 
around a constant value, the external increase must be exactly compensated by the 
average dissipation. This gives an exact formula involving the average earthquake size 
(s) defined as the mean number of discharges per earth quake, 

{l-na)z;^{s) = N{l-z^). (3) 

Notice that the product of averages on the l.h.s. does not result from a factorization 
approximation but from the definition of Zunst ? and is exact. Therefore, this equation 
is correct even if the above mentioned simplifying assumptions are not true, and holds 
thus also in the fixed neighbor version of the model. Since the left hand side of this 
equation remains finite for iV — >• cx) (as long as a < 1/n), we see that 1 — z^ oc 1/N. 

On the other hand, since the force increase between earth quakes is assumed to be 
with velocity v = 1, the average time between two earth quakes is given by 1 — z^. On 
a 'macroscopic' time scale where we neglect the duration of earth quakes compared to 
the inter-quake times (this assumption is inherent in the model), the toppling rate is 
thus given by 

(s) ^ 1 
N{l-z^) (1-na)^;;^' 

This tells us how frequently each site discharges per time unit. The rate to be hit 
according to eq.(0) is then given by na. 

Let P{z) be the probability density for a given site to have a force value z (from 
now on, we shall consider only = oo). Obviously, -P(l) is the rate with which new 
earth quakes are initiated, while -P(O) = a is the rate with which new force-free sites 
are created by discharges. Therefore, 

(.) = P(0)/P(1) . (5) 

Similarly, Pj{z) denotes the joint probability density that a site has a value z and 
was hit exactly j-times since its last discharge. Since we consider only N = oo and 
have argued that global avalanches are negligible in this limit, P{z) and Pj{z) do not 
fluctuate with time. Obviously, we have 

m 

Piz) = Y.P,{z) ze[0,l]. (6) 

j=Q 

Because a hit increases z at least by an amount a, each Pj{z) vanishes exactly for 
z < ja. Therefore the upper limit m in the above sum is given by the largest integer 
for which ma < 1. For later use we define the integrated distribution as 

V{z) = C P{z')dz' . (7) 

J z 

To obtain Po{z) we notice that the probability to be hit exactly k times during a 
time interval z, when the rate is na, is given by the Poisson distribution 

Mz) = ^(n^^)' e""'^' . (8) 
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This leads to 

Poiz) = anoiz) = ae-'^'^^ (9) 

The other Pj{z) depend on the distribution of the amount Az which a site receives 
when it gets hit by a discharge. This in turn depends on the distribution of forces 
of unstable sites at the moment of their discharge. We denote the density of this 
distribution by C (z) . It is related to z^ by 

/oo 
z C{z)dz. (10) 

(Here and in the following, integrals over functions with (5-peaks at the integration 
limits are understood as containing all contributions from these peaks, Jl^f{x)dx = 
lim^^Q JI^^^ f{x)dx.) The first site of any earthquake discharges exactly with z = 1. 
This gives a delta contribution to C{z), with relative weight l/(s). We can therefore 
make the ansatz 

C{z) = 1^5{z-1) + C{z). (11) 

The second term corresponds to all subsequent discharges. About the function C{z) 
we know that it has to vanish for all z outside the interval [1, 1/(1 — a)]. The upper 
limit would be reached if an infinite earthquake contained a series of successive hits 
onto sites with maximum value z = 1. This upper limit could be surpassed only if a 
site were hit simultaneously by two discharges, but the chance for this is zero on an 
infinite lattice. The amount of force Az that a discharging site drops onto each of the 
n random neighbors is then distributed according to 

Qi{Az) = a~^C{Az/a) , supp Qi = [a, a/(l - a)] . (12) 

Similar to eq. (pA]) we can write Qi{Az) as 

Qi{Az) = ^5{Az-a)+Qi{Az). (13) 

The convolution integrals 

QkiAz) = Qk-ii^z - Az')Qi{Az')dAz', k>2 (14) 

J a 

then give us the probability densities for the total increase of force when a site was hit 
exactly k times. Note that Qk{Az) vanishes for Az outside the interval [ka, ka/{l—a)]. 
We see finally that every Pj{z) has to obey 



PoiO) r 7ij{z- Az)Qj{Az)dAz 

Jja 

r {na{z - Az)y e-'"'^'-^'^Qj{Az)dAz. (15) 

J ja 



a 



Let us now come back to the function C{z) which describes the distribution of 
sites which have become unstable by being hit by a discharge. It is obtained by the 
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convolution of the distribution of (stable) sites before they are hit, with the distribution 
of amounts received during the destabilizing discharge, 



C{z) =nQ{z-l) 



P{z - Az)Qi{Az)dAz 



(16) 



Here Q{x) is the Heaviside step function. It takes into account that C{z) is supposed 
to describe only those sites which actually do become unstable. The factor n takes into 
account that each discharge event — the probability density of which is given by the 
integrand — gives rise to n potentially unstable sites. The integration limits are given 
by the support of Qi, see eq.(|T^). 

Let us now consider the integral of C{z) over all z. Interchanging the integrations 
we obtain 

/ C{z)dz = n dAzQi{Az) P{z - Az) dz . (17) 

Jl Ja Jl 

The inner integral on the rhs. is just V{1 ~ Az) (see eq.(|^)), while the left hand side 
is equal to 1 — l/(s) due to eq. (pJ]) . Rearranging terms we arrive thus at a second 
equation for (s), 

(s) = 1-n r^''^ V{1- z)Qi{z)dz . (18) 

J a 

We claim that the above equations are complete in the sense that they fix the 
solution uniquely, for each a < 1/n. To show this, and to provide also a practical 
method to solve them numerically for not too large a, we give a recursion scheme 
which converges to the solution as the iteration level r tends to infinity, at least for 
sufficiently small a. For larger values of a the recursion might not be practical, but 
the set of equations should still fix the solution by continuity. Notice that different 
recursion schemes are in principle possible where the order of replacements is changed 
in various places. 

To start the recursion, we select a desired accuracy rj and choose the initial distri- 
bution Qi{zY^^ in some arbitrary way. It need not even be normalized. For small a, a 



good choice is Qi{z)^'^^ constant. For a ~ 1/n we can also take Qi{z)^^^ 
In the recursive step from r — 1 to r we do the following: 

(1) (re-) normalize Qi{z) = const x Q^^ with const = [J ^\z)dz]^^ 

(2) compute a from eqs.(|D, (|10D, and (|l|), 

-1 



6{z-l/n) 



a 



a 



'1 — na) 



a/{l-a) 



Qi{z) z dz 



(19) 



(3) compute Qk{,z) for > 1 by means of eq.(|T^ 



(4) compute Po{z) = ae compute Pk{z) for /c > by means of eq.(|T5|), and obtain 
P{z) as Y.kPk{z)- 

(5) compute the new (s) from eq.(||); 

(6) compute the new Qi{z) from eqs.(|T2|), (pA]), and (|16|), 

ra/(l-a) 



a) + - Qiz 
a 



a 



P{zla-OQ,{C)dC 



(20) 
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(7) if 0" or (s) have changed by a fraction larger than t], then goto (1); 

(8) verify that the normahzation constant in step (1) is unity within some acceptable 
error, and that (s) satisfies eq. ([T8D . 

We have not shown formally that this iteration converges always, but we have done 
extensive numerical investigations. The scheme converges very fast and stably for small 
a, but convergence is slowed down when a 1/n. For this reason we had problems 
to obtain solutions for a > 0.24, although the recursions shows no sign for divergence 
even in such extreme cases. Also, numerical errors in the integration routines tend 
to accumulate for a 1/n, rendering in particular the estimate of (s) problematic. 
Since the integrands are not analytic functions, it does not make sense to use very 
sophisticated integration routines. We used the extended trapezoidal rule with up to 
10'' points. 

12 I 1 1 1 1 1 

10 - 
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Figure 1: Probability density P{z) against z. The continuous lines are the predictions from 
the theory and the points show the results obtained from simulations. They fall perfectly on 
top of each other on the scale of this figure. The four curves are for a = 0.20, 0.21, 0.22, and 
0.23, in order of increasing sharpness of the peaks. 



Results for n = 4 are shown in fig. |l], where we also compare with straightforward 
simulations of eqs. (|1|),(|). For the latter we typically used = lO'' to 8 x 10^ and 
discarded transients of up to 2 x 10^ iterations. No difference between theory and simu- 
lation is detectable in fig. |l|. This shows that the numerical integration was sufficiently 
accurate, the iterations had converged, N was sufficiently large to have negligible finite 
size corrections, and the discarded transients were sufficiently long. Qualitatively, fig. |] 
is similar to fig.l in [IT^, but the first peak in that paper seems much too high. It is 



not clear whether this results from the be used in that paper or from transients. For 
a = 0.23, we find -P(O) ~ 11.5 both from simulations and from the analytic solution, 
while P(0) ^ 33 is quoted in [0;. 

In fig. II we show a and (s) as functions of e = 1 — 4a. We see that a ~ 1/e, 
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Figure 2: Log-log plot of a and (s) against e = 1 — Aa. Continuous lines are from theory, 
points from simulations. 



while (s) diverges much faster when a — >■ 1/4. Finally, in fig. ^ we show Qi{z). This 
shows a very interesting qualitative change as a approaches 1/4. For a < 0.23, Qi{z) 
is centered at 2; < 1/4. Its center moves to the right as a increases, reaching a value 
slightly larger than 1/4 for a ~ 0.233. After that, its center moves very little, and it 
just shrinks slowly to a 5-function centered at 1/4. 

The most important result is that we see no hint of any singularity for o; < 1/4, as 



predicted in [|I^ , and we also see no mechanism which could lead to such a singularity. 
Indeed, we can prove rigorously that (s) < 00 for all a < 1/n. This follows simply 
from the fact that a < 1/(1 - na) due to eq.(4), and (s) = a/P{l) < ct/Po(1) = e"'' 
due to eqs.(5) and (9). Conversely, this argument shows that P{1) must tend to for 
e — > 0, as also shown by the numerics. 

According to |1^, a singularity with (s) — > 00 should occur for ?T, = 4ata = 2/9 = 



0.222 .... We beheve that this is due to unjustified assumptions made in [T2|. Another 



important result is that P{z) is finite and non-zero at z = 1 and at 2; = 0. This shows 
that 'global' earth quakes have indeed no effect, as they would lead to a depletion at 
P{z) at 2 = 1 or an infinity at 2 = due to eq.(^). 



4 The Limit a ^ 1/n 

Figure 1 suggests that P{z) tends to a sum of four delta peaks at multiples of a, for 
a 1/4. More generally, we expect P{z) to tend towards a sum of n delta peaks 
at z = k/n, k = 0, . . .n — 1 (this is reminiscent of a generalized sandpile model with 
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Figure 3: Qi{Az) against Az, for the same values of a as in fig. |l[ and in addition for 
a = 0.233 and a = 0.236. For clarity only theoretical predictions are shown. The cusps are 
at Az = a + ct^ and correspond to the maximal Az transfered by first generation descendents 
of the avalanche seed. 



real- valued heights by Zhang [Q). We shall see that this is indeed a valid solution 
after proper rescalings of a and (s). 
Formally we introduce 

e = 1 - an , (21) 

and consider the limit e — >■ 0. We shall argue that a self consistent solution for P{z) 
in this limit is a sum of delta peaks. If this is true, only sites which have already 
z ^ {n — l)/n will become unstable by receiving an extra Az ^ 1/n, and hence 
-^unst 1 for e ^ 0. From this we see on the one hand that a = 1/e to leading 
order in 1/e, which in turn gives aTikiz) n^^6{z) for each value of k. On the 
other hand it gives C{z) = S{z — 1) and Qi{Az) = 6{Az — 1/n). The latter implies 
Qk{Az) = 6{Az — k/n) for any k > 2, which finally gives 



n-l 



P{z) = -J26{z-j/n), 



(22) 



n 



j=0 



i.e. our initial assumption was self consistent. 

In spite of the simplicity of this solution, we should be careful in interpreting it, as 
several limits are involved. The easiest way is to take first the infinite volume limit, 
and then e — 0. If we want to take the limit e ^ first, we have to use absorbing sites 
which mimic absorbing bc's, but it is not a priori clear how their number should scale 
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with (for a related problem in a mean field version of the abelian sandpile model, 
see 0,111). 

While the behavior exactly at a = l/ra is thus well understood, we were not able to 
find an analytic solution for finite e. But we can give approximate solutions for small 
e, and predict the behavior of (s) for e — >■ 0. For small but finite e, we approximate 
each Qj{z) by a delta function at z = j/n. Then Pj{z) is roughly given by 

P,{z) ^ e{z-j/n) ^[{z-j/n)aye-^^-^/"^^ . (23) 
From this, eq.(|D, and the fact that -P(O) = cr, we get 

(.) = a/P{l) ^ a/P„_i(l) ^ ^^^e'^ ^ {n - l)\e--'e'/^ . (24) 

There are substantial corrections to this, mainly from the contribution of -Pn(l) to -P(l), 
which are hard to estimate. Thus, the actual values of (s) are smaller than given by 
eg. (p^ , but eq.(p^) gives the correct trend. In particular, it explains why (s) diverges 
extremely fast for a — >■ 1/n, making simulations in this limit very difficult. 



5 Earth Quake Statistics 

In the previous sections we have only studied the force distribution and the average 
earth quake size. In order to discuss the distribution of earth quake sizes and durations, 
we need some more definitions and some basic results from the theory of branching 



processes as found e.g. in [|T7 



For any integer i > 1 we define Pi as the probability that a site becomes unstable 
if it is hit by a discharge event during the {i — l)-st generation of an earth quake, 
and will therefore discharge itself in the i-th generation. For the first generation, pi 
is the probability that a site which receives Az = a gets unstable. Thus, simply pi = 
V{1 — a). For general a, the other probabilities Pi depend on the height distribution 
of the unstable sites in the previous generation. But for e — >■ all discharging sites 
have z = 1, and pi is simply the probability that the hit site is in the (n — l)-st peak. 
Pi = p = 1/n for alH > 0. 

In the following we shall therefore discuss only the case e = 0, deferring the general 
case to the end of this section. 

We assume thus that Pi = p = 1/n for all i > 0. The probability that an unstable 
site creates / unstable offsprings in the next generation is then given by 



wi=^^jil-pr-V (25) 

with the generating function 

n 

g{u) = Y,u^wi = (l-p + up)''. (26) 
z=o 
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With the use of g{u) it becomes very easy to calculate recursively the distributions 
of the size and the duration of the earthquakes. In order to do this we introduce the 
iterates 

9o{u) = u, gi{u) = g{u), 
gm+iiu) = g[gmiu)], m = l,2,...,. (27) 

The probability that the earthquake stops in the t-th time step is given by JT^ 

Vt = gt{0) - gt-i{0) with t=l,2,.... (28) 




Figure 4: The integrated distribution Pj. Again, lines show the theoretical predictions and 
points are from simulations. Increasing from left to right, a takes the same values as in fig. 1. 



The dashed line shows eq. (32). 



The integrated distribution 



(29) 



t'=t 



denotes the probability that an avalanche lasts for > t time steps. Eq. (^91) can be used 
directly to calculate Pt for small t, whereas for large t we use the asymptotic behavior. 
The chain of identities 

Pt+i = l-gt{0) 

= l-g{gt-im 
= l-g{l-P,) 
= l-(l-pPi)" 



(30) 
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leads to 

^ P,+i - = - hp^P,^ + o{P!) (31) 



with the solution 



Pt -t-^ ■ 32 

n — 1 



This is a special case of the general theorem |]T 

2 



(33) 



for a critical branching process. 

The next quantity of interest is the size distribution of the earthquakes. With T>s we 
denote the probability that the size of an avalanche is exactly s. While Vi = (1 — 1/n)" 
is obvious, the calculation of Vg for s > 1 proceeds as in |T8[. We first denote by a^^^ 
the k-th Taylor coefficient of [(^(m)]*, i.e. [(7(^)1^ = a^''' + ai'^ u + a^2^ u'^ + . ... A theorem 
due to Dwass p2| tells us then that 

V, = ^ai'l,, s>l. (34) 

In the present case, we have 

4^^ = (7)(i-pr"V, (35) 



leading to 



Vs = -(''"\l-p)^-~'>^y-\ (36) 



The local limit theorem of Moivre-Laplace states that in the limit s ^ oo the distri- 
bution Vs with p = 1/n tends to 

^ ^ ^ s-'/\ (37) 

This means that the system gets critical for a = 1/n, and the critical exponents take 
the same values as for mean-field percolation |2l[ . 



In the subcritical phase the arguments are more tedious. Let us denote by Ci{z) 
the joint probability distribution that a discharge happens during the i-th generation 
of an earth quake, and that the discharging site has force z. It is related to C{z) and 
to Pi by 

C{z) = -^S{z - 1) + c,{z) + C2{z) + .... (38) 

and 

„ . ii;;;:'-'-)"- . (39, 

n Ji Ci_i{z) dz 
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In s 

Figure 5: The integrated distribution Dg = X]^=s ^s' with the same a-values as in fig. ^. 
The dashed hne shows the scaling law = 2 g-^/^. 

This is easily checked by noting that it is compatible with 

(s) = 1 + npi + n^piP2 + n^PiP2P3 + • • • • (40) 
The functions Ci{z) satisfy a recursion relation similar to eg. ([T6|) , 

rl/{l-a) 

Ci{z) = n Q{z - 1) P{z-az')ci^i{z')dz (41) 

with Cf){z) = 6{z — 1)/ (s). Again we were not able to solve this analytically. But given 
a numerical estimate of P{z) we can solve it numerically for Ci{z), i = 1,2, . . ., from 
which we obtain pi by integration. Again this was done only for n = 4. For each 
considered value of a we found that Pi increases monotonically with i and converges 
very quickly to a constant value < 1/4. This is easy to understand. The increase is 
due to the fact that the first discharging site has z = 1, while all subsequent ones have 
z slightly larger than 1. The fact that pi < 1/4 reflects the fact that we are dealing 
with a subcritical branching process. 

In contrast to the critical case we now have a time dependent (non-autonomous) 
branching process, i.e. the generating function g{u) depends on the generation. The 
mathematical treatment becomes now more tedious. For the further comparison be- 
tween simulations and theory we therefore used the theoretically obtained pi to sim- 
ulate a branching process, and compared the results with direct simulations of the 
OFC model. Figures § and ^ show that the agreement is essentially perfect, except for 
large s and t, and for a = 0.23. The discrepancies seen there arise from the numerical 
problems mentioned in sec. 3. 
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6 The Feder-Feder Model 



Up to now we have considered only the OFC model, but our methods are much more 
general. To illustrate this, we shall discuss in this section the random-neighbor version 
of a model introduced by Feder and Feder (FF) [Q. The FF model is identical to 
OFC model, with the only exception that eq.(|D is replaced by 

^ik ^ + a . k = l,...,n. (42) 

This means that a site hit by a discharge always receives a fixed amount a regardless 
of the z-value of the discharging site. This leads to a significant simplification of the 
equations of motion. To derive them, we have first of all to notice that the toppling 
rate is now given by 

..(fl^^^—. (43) 

N{l-Zm) Zunst-na 

Since Zunst > 1, it is not a priori clear whether a diverges in the limit a 1/n. But 
we will show that this is indeed the case. 

The main simplification arises from the fact that Qi{Az) uncouples from C{z) and 
is given by 

Qi{Az) = 5{Az - a) . (44) 
From this one obtains immediately 

QkiAz) = SiAz - ka) , (45) 

see eq.(|l^, and 

m /T 

Pi^) = E^.(^) = E -ina{z - ja)ye{z ~ ja)e-'^''^'-^''\ (46) 

where m is again the largest integer < 1/a. To obtain a as a function of a, we use the 
normalization condition P{z)dz = 1 which gives 

™ 1 /•{l—aj)na 

n = y- dxx^ e"^ . (47) 

For a close to 1/n we have m = n, and this condition can be rewritten as 

dx e"" = V - / dx x' e"" , (48) 

^-^0 j! i{l-aj)na 

where we have used again e = 1 — na. From this it is easily seen that a — oo for 
e — > 0. Otherwise, the left hand side would tend to zero in this limit, while the r.h.s. 
would remain non-zero. But if a diverges, the r.h.s. is dominated by the term with 
j = n — 1. Keeping only dominant terms we arrive at 

a ^ (n + 1) log(l/e) . (49) 
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The mean avalanche size can again be calculated from eq.(^, 



(s) = a/P(l) 



^7rj (l - ja) 
i=o 



(50) 



Obviously (s) is finite for all a < 1/n and diverges for a ^ 1/n. In this limit the sum 
is dominated by the term with j = n, giving 



(s) ~ n\ {naeY 



(51) 



up to constant and logarithmic factors in e which could easily be computed. 

Thus cr and (s) both diverge much slower than in the OFC model. This reflects the 
increased dissipation in the FF model. Exact values of a and (s) obtained numerically 
from eqs.(^) and ( ^OD are shown in fig.6. For small e one finds good agreement with 
the asymptotic predictions. Avalanche dynamics can be treated exactly as in the OFC 
model. 
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Figure 6: Log-log plot of a and (s) against e = 1 — 4q for the Feder-Feder model. 



7 Conclusion 

Our results show clearly that there are neither scaling nor phase transitions in the 
random neighbor version of the dissipative Olami-Feder-Christensen earthquake model. 
Scaling is observed only in the conservative limit, in which case one has a critical 



branching process. This is in direct contradiction to claims in ||1^. The latter was 
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based on approximate random neighbor equations, while the present work is based 
on the exact equations. These equations were solved numerically, giving excellent 
agreement with direct simulations of the model. 

The most surprising result was the very fast increase of the average earth quake size 
as one approaches the conservative limit. Obviously this is a consequence of the non- 
locality of the interaction, since this implies that one can have extremely large earth 
quakes without having large effects locally. Nevertheless, avalanche size distributions 
decay exponentially for any non-zero dissipation. 

Our findings support the view that scaling in the OFC model with inhomoge- 
neous boundary conditions is due to a subtle interplay between partial synchronization 
and desynchronization. The inhomogeneity of the be drives the synchronization in the 
bulk, building up large coherent patches, but occasionally the driving is too strong and 
the synchronization breaks down. Explicit observations of these patterns ||^, § support 
this view. In a random neighbor model such structures cannot build up, of course, and 
the mechanism driving the system into a self-organized critical state is absent. 

While we concentrated here on the OFC model, we showed that our methods can 
also be applied in other related models. In particular, we studied the Feder-Feder model 
in some detail. We showed that it also has no phase transition in the dissipative regime, 
and that the toppling rate and the mean avalanche size diverge in the conservative limit. 
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